{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "利用GDAL实现波段运算和阈值分割"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from osgeo import gdal \n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_img=\n",
    "dset=gdal.Open(path_img)\n",
    "img_array=dset.ReadAsArray()\n",
    "print(img_array.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "band_green=dset.GetRasterBand(2).ReadAsArray().astype(np.float16)\n",
    "band_swir=dset.GetRasterBand(5).ReadAsArray().astype(np.float16)\n",
    "mndwi=(band_green-band_swir)/(band_green+band_swir+0.00000000001)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10,3))\n",
    "plt.subplot(1,3,1)\n",
    "plt.imshow(band_green,vmax=3000,vmin=0,cmap='gray')\n",
    "plt.subplot(1,3,2)\n",
    "plt.imshow(band_swir,vmax=3000,vmin=0,cmap='gray')\n",
    "plt.subplot(1,3,3)\n",
    "plt.imshow(mndwi,vmax=1,vmin=-1,cmap='gray')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "id_valid=np.where(band_green>0)\n",
    "\n",
    "n,bins,patches=plt.hist(mndwi,bins=150)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "thre=0.25\n",
    "mndwi_thre=np.where(mndwi>thre,1,0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path_out=\n",
    "driver=gdal.GetDriverByName('GTiff')\n",
    "out_dset=driver.Create(path_out,xsize=dset.RasterXSize,ysize=dset.RasterYSize,\n",
    "                   bands=1,eType=gdal.GDT_Int16)\n",
    "out_dset.SetGeoTransform(dset.GetGeoTransform())\n",
    "out_dset.SetProjection(dset.GetProjection())\n",
    "out_band=out_dset.GetRasterBand(1)\n",
    "out_band.WriteArray(mndwi_thre)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.18"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
